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Abstract. Correlations are employed in modern physics to explain microscopic 
and macroscopic phenomena, like the fractional quantum Hall effect and the 
Mott insulator state in high temperature superconductors and ultracold atoms. 
Simultaneously probed neurons in the intact brain reveal correlations between 
their activity and theoretical work illuminates the importance of correlation for 
information processing and for macroscopic measures of neural activity, like the 
EEC Nevertheless networks of spiking neurons differ from most physical systems: 
The interaction between elements is directed, time delayed, is mediated by short 
pulses, and each neuron receives events from thousands of neurons and sends events 
to thousands of others. Even in the stationary state the network does not reach 
equilibrium in the sense of detailed balance. Here we develop a quantitative theory 
of pairwise correlations in finite sized random networks of spiking neurons. We 
show why the intuitive mean field description fails, how single action potentials 
reverberate in the network causing an apparent lag of inhibition with respect to 
excitation, and how global collective oscillations arise. 
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1. Introduction 

Correlations are an established feature of neural activity [1] and evidence increases 
that they are essential for information processing in the brain [2]. The temporal 
relationship between the activity of pairs of neurons is described by correlation 
functions. Their shape has early been related to the direct coupling between neurons 
and to the common input shared by pairs of neurons. Correlations affect the 
representation of information as they may limit the signal-to-noise ratio of population 
rate signals [3]. On the other hand, they have been shown to increase the amount 
of information available to unbiased observers [4]. Furthermore, synchronous neural 
activity is thought to bind elementary representations into more complex objects [5] 
and experimental evidence for such a correlation code is provided by task related 
modulation of synchrony in primary visual cortex [6] and in motor cortex [7]. 

The small magnitude [8] of pairwise correlations in the asynchronous irregular 
state [9] of cortex has recently been related to the balance between excitation and 
inhibition in local networks [10, 11] and inhibitory feedback was identified as a 
general mechanism of decorrelation [12]. However, a quantitative theory explaining 
the temporal shape of correlation functions in recurrently impulse coupled networks 
of excitatory and inhibitory cells remained elusive. 

Assuming random connectivity with identical numbers and strengths of 
incoming synapses per neuron, as illustrated in Figure 1, suggests by mean field 
arguments [13, 14] that the resulting activity of two arbitrarily selected neurons 
and hence the power spectra of activities averaged over excitatory or inhibitory 
neurons should be the same. Direct simulations, however, exhibit different power 
spectra for both sub-populations [15]. A similar argument holds for the covariance Cfj 
between the two neurons: If the covariance c between any pair of inputs is known, the 
covariance between their outgoing activity cg is fully determined. By self-consistency, 
as both neurons belong to the same recurrent network, one concludes that cg = c. 
In particular the covariance averaged over excitatory pairs should be identical to the 
corresponding average over inhibitory pairs, which is in contrast to direct simulation 
(Figure lb). In this work, we elucidate why this mean field argument for covariances 
fails and derive a self-consistency equation for pairwise covariances in recurrent 
random networks which explains the differences in the power spectra and covariances. 

Theories for pairwise covariances have been derived for binary neuron models 
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Figure 1. Self-consistency argument fails for covariances in a homogeneous 
recurrent random network, (a) Each neuron (black circles) receives input from 
the same number of randomly chosen excitatory (e) and inhibitory (i) neurons in 
the network, so the input statistics of all neurons is the same. The covariance c 
within the network determines the covariance between the neuron inputs to a pair 
of neurons and hence the covariance eg of their outputs. Self-consistency seems 
to require eg — c = c^c = ca. (b) Covariance functions averaged over pairs of 
excitatory (cqg) and over pairs of inhibitory (ca) integrate- and- fire model neurons 
are different in a direct simulation. The other parameters are given in Appendix E. 

[11, 16] and for excitatory stochastic point process models [17]. However, the lack 
of either inhibition or delayed pulsed interaction limits the explanatory power of 
these models. A theory for networks of leaky integrate-and-fire model neurons [18] 
is required, because this model has been shown to well approximate the properties 
of mammalian pyramidal neurons [19] and novel experimental techniques allow to 
reliably assess the temporal structure of correlations in cortex [20]. Moreover, the 
relative timing of action potentials is the basis for models of synaptic plasticity 
[21], underlying learning in biological neural networks. Analytical methods to treat 
population fluctuations in spiking networks are well advanced [22] and efficient hybrid 
analytical-numerical schemes exist to describe pairwise covariances [23]. Here we 
present an analytically solvable theory of pairwise covariances in random networks 
of spiking leaky integrate-and-fire model neurons with delayed pulsed interaction in 
the asynchronous irregular regime. 

2. Results 

We consider recurrent random networks of excitatory and -"yN inhibitory leaky 
integrate-and-fire model neurons receiving pulsed input (spikes) from other neurons 
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in the network. Each neuron has K = pN incoming excitatory synapses 
independently and randomly drawn from the pool of excitatory neurons, and 
'jK = '-fpN inhibitory synapses (homogeneous Erdos-Renyi random network with 
fixed in-degree). An impulse at time t arrives at the target neuron after the synaptic 
delay d and elicits a synaptic current Jj that decays with time constant and causes 
a response in the membrane potential Vi (with time constant Tm) proportional to the 
synaptic efficacy J (excitatory) or —gj (inhibitory), respectively. The coupled set 
of differential equations governing the subthreshold dynamics of a single neuron i is 
[24] 

dVi 

Tm-^ = -Vi + h{t) 

dl ^ 

Ts-^ = - li + JijSj{t - d), (1) 

where the membrane resistance was absorbed into Jij. If Vi reaches the threshold Ve 
at time point t]. the neuron emits an action potential and the membrane potential is 
reset to Vr, where it is clamped for the refractory time r^. The spiking activity 
of neuron i is described by this sequence of action potentials, the spike train 

The activity of a given neuron i depends on the history of the other neurons' 
activities s(t) = {si(t), . . . , si\[(t))'^ in the network. The time averaged covariance 
matrix expresses these interrelations and is defined as c(r) = ^s(t + r)s-^(t))^ — rr'^, 
where r = (s)^ is the vector of time averaged firing rates. The diagonal contains the 
autocovariance functions (diagonal matrix a(r)) which are dominated by a 5-peak 
at zero time lag and for r 7^ exhibit a continuous shape mostly determined by 
refractoriness, the inability of the neuron to fire spikes in short succession due to the 
voltage reset, as shown in Figure 2d. The off-diagonal elements contain the cross 
covariance functions c(r) that originate from interactions. We therefore decompose 
the covariance matrix into c = a + c. A basic property of covariance matrices is the 
symmetry c(r) = c^{—t), so we only need to consider r > and obtain the solution 
for r < by symmetry. Each spike at time t' can infiuence the other neurons at time 
t > t' . Formally we express this infiuence of the history up to time t on a particular 
neuron i in the network as Si{t, {s{t')\t' < t}), which is a functional of all spike times 
until t. In the asynchronous state of the network [9] and for small synaptic amplitudes 
Jij a single spike of neuron j causes only a small perturbation of Sj. The other inputs 
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to neuron i effectively act as additional noise. We may therefore perform a linear 
approximation of j's direct influence on neuron i and average over the realizations 
of the remaining inputs s\sj, as illustrated in Figure 2a. This linearization and the 
superposition principle lead to the convolution equation 

{s^{t)\s,Us,=r,+ I K,{t,t'){s,{t')-r,)dt' (2) 

J —oo 

= ri + [hij*{sj-rj)]{t), (3) 
where we define as the linear response kernel hij (t, t') = ( ss^W) / functional 

\ ' / s\Sj 

derivative of Si(t) with respect to Sj(t'), formally defined in Appendix A. This kernel 
quantifies the effect of a single spike at time t' of neuron j on the expected density of 
spikes of neuron i at time t by a direct synaptic connection from neuron j to neuron 
i. This density vanishes for t < t' due to causality. For the stationary network state 
studied here, the kernel further only depends on the time difference r = t — t' . In a 
consistent linear approximation there are no higher order terms, so the effects of two 
inputs Sj and Sk superimpose linearly. Using this linear expansion and the definition 
of the covariance matrix, the off-diagonal elements of the covariance matrix fulfill a 
linear convolution equation 

c(r) = [h* (a + c)](r) for r > 0. (4) 

The equation is by construction valid for r > and needs to be solved simultaneously 
obeying the symmetry condition c{t) = c(— t)^. Up to linear order in the interaction, 
equation (4) is the correct replacement of the intuitive self-consistency argument 
sketched in Figure 1. 

In order to relate the kernel hij to the leaky integrate-and-fire model, we 
employ earlier results based on Fokker-Planck theory [24]. For small synaptic 
amplitudes J <^Ve — Vr and weak pairwise covariances the summed synaptic input 
Tm'^f^i JijSj{t — d) can be approximated as a Gaussian white noise with mean 
/Wj = Tm J2j Ji^j ^'^d variance af = t„i J^fj neurons firing with rates and 
Poisson statistics. For short synaptic time constants ^ the stationary firing 
rate ri{fj,i,af) (A.l) depends on these two moments [24]. In a homogeneous random 
recurrent network the input to each neuron is statistically the same, so the stationary 
rate = r is identical for all neurons. It is determined by the self-consistent solution 
of r(/ij, af), taking into account the dependence of /ij and af on the rate r itself [9]. 

The integral Wij = hij (t) dt of the response kernel is equivalent to the DC 
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Figure 2. Mapping the integrate- and- fire dynamics to a linear coupling kernel, (a) 
The kernel /i^ determines the transient effect of an incoming impulse at time point 
tj (black arrow) on the density Si{t) of outgoing action potentials, averaged over 
realizations of the stochastic activity of the remaining inputs (indicated as red and 
blue triangles), (b) Response kernel (5) (green) compared to direct simulation for 
an impulse of amplitude J = 1 mV (black dots) and J = — 1 mV (gray dots). Time 
constant ~ 4.07 ms determined by a least squares fit to a single exponential. The 
background activity causes a mean fii = 15 mV and fluctuations ai = 10 mV. (c) 
Linear and quadratic dependence (A. 3) of the integral response Wij on Jij (dark 
gray curve) and linear term alone (light gray line), (d) Autocovariancc function of 
the spike train with a 6 peak at i = and covariancc trough due to refractoriness. 



susceptibility Wij = [25], leading to an approximation of second order in the 
synaptic amplitude Wij = a Jij + /3Jj^-. The first order term originates from the 
dependence of /ij on Vj, the second order term stems from the dependence of af on 
rj (A.3). 

Figure 2b shows the deflection of the firing rate from baseline caused by 
an impulse in the input averaged over many repetitions. For sufficiently strong 
fluctuations (Tj, the quadratic term in Jij is negligible, as seen from Figure 2c. The 
kernel shows exponential relaxation with an effective time constant Tg, that depends 
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on the working point (/ij, (Tj) and the parameters of the neuron, which is obtained by 
a least squares fit of a single exponential to the simulated response in Figure 2b for 
one particular amplitude Jij. We therefore approximate the response hij{t) as 

WiMt) = e{t - d) ^e"^, (5) 

Te 

where d is the synaptic delay and 6 the Heaviside step function. 

In experiments covariance functions are typically averaged over statistically 
equivalent pairs of neurons. Such averages are important, because they determine 
the behavior of the network on the macroscopic scale of populations of neurons. We 
therefore aim at a corresponding effective theory. Assuming identical dynamics (1), 
all neurons have to good approximation the same autocovariance function a{t) and 
response kernel h{t). So each incoming excitatory impulse causes a response wh{t), 
an inhibitory impulse —gw h{t). We define the covariance function averaged over all 
pairs of excitatory neurons as Ccc(t) = -^J^ijeeij^j^ij^'^)^ (setting A^(A^ — 1) ^ A^^ 
for N ^ 1), where S denotes the set of all excitatory neurons. The pairings Cei,Cie, 
and Cii are defined analogously. Inserting equation (4) into the average Ccc{t), 
the first term proportional to the autocovariance a{t) only contributes if neuron 
j projects to neuron i. For fixed i, there are K such indices j, so the first term 
yields Y.i,j(i£,i^j ^ij * dj = NKw h * a. The second sum Xli jef * Ckj can be 

decomposed into a sum over all intermediate excitatory neurons k E £ and over all 
inhibitory neurons k eX projecting to neuron i. Replacing the individual covariances 
by their population average, Cee and Cie, respectively, and considering the number of 
connections and their amplitude we obtain wNKh * (cee — '^gcic). Similar relations 
hold for the remaining three averages, so we arrive at a two-by-two convolution 
matrix equation for the pairwise averaged covariances for r > 



c(r) = [/i * Mc] (r) + Q [/i * a] (r) (6) 

Kw 



with m =Kw\ ^ 1 > Q ^ 




and c{t) 



Cicir) Cii(r) 



The convolution equation only holds for positive time lags r. For negative 
time lags it is determined by the symmetry c(-r) = c^(r). The solution 
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of this equation can be obtained by an extension of the method used in [17] 
employing Wiener-Hopf theory [26] to the cross spectrum C{u) = c(t)e"*'^* cit 
in frequency domain, as shown in Appendix B (here capital letters denote the 
Fourier transform of the respective lower case letters). With the definition 
of the propagator P{u)) = (1 — Mi7(a;))~^ the cross spectrum takes the form 
C{uj) = P(w) [D4uj)q + Z)+(-cj)Q^ - A{uj)\H{uj)\^QM^] P^{-uj), where we split 
the term H{uj) A{u) = D^{uj) + D_{uj) so that (i+(r) and (i-(r) vanish for times 
r < and r > 0, respectively. The non- averaged cross spectrum can be recovered 
as a special case setting Q = M = W, because the convolution equations (4) and 
(6) have the same structure and symmetries. If all eigenvalues of Mif(u;) have an 
absolute value smaller unity, the propagator P{uj) can be expanded into a geometric 
series in which the n-th term contains only interactions via n steps in the connectivity 
graph. This expansion has been used to obtain the contribution of different motifs 
to the integral of covariance functions [27] and their temporal shape [28]. 

In the following, we neglect the continuous part of the autocovariance function, 
setting a{t) = r S{t), because the 6-peak is typically dominant. For an arbitrary 
causal kernel h it follows that D^{u) = r H{u), so the cross spectrum takes the form 

C{uj) = r^ \ ^^^^^ ^ ^^"^ 

and L =Kw{l — 'yg). 

The limit w — corresponds to the time integral of the cross covariance function, 
approximating the count covariance for long time bins [2, 8]. With A{0) = r, the 
integral correlation coefficient averaged over neuron pairs fulfills the equation 



C(0) _ Kw 1 ( 2 l-g\ 

{Kw^l+g^ ( 1 l\ 
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Figure 3. Shape invariance of covariances with network scale. Synapses scaled 
as J cx 1/A^ with network size N to conserve the population feedback L = const, 
(a) Integral correlation coefficient averaged over different pairs of neurons and 
theory (8) confirming the (x 1/A^ dependence, (b) Rescaled covariance functions 
Ncei averaged over excitatory-inhibitory pairs of neurons for different network sizes 
(color coded) and theory (15) (black). 



which has previously been derived from a noise-driven hnear rate dynamics [12]. 
The quantity L plays a key role here: It determines the feedback magnitude of 
in-phase fluctuations of the excitatory and inhibitory population. Stability of the 
average firing rate requires this feedback to be sufficiently small [9], i.e. L < 1, 
indicated by the pole in equation (8) at L = 1. Typically cortical networks are in 
the balanced regime [9, 13, 14], i.e. L < 0. For such inhibition dominated networks, 
the denominator in equation (8) is larger than unity, indicating a suppression of 
covariances [12]. As shown in Figure 3a, the prediction (8) agrees well with the results 
of simulations of leaky integrate-and-fire networks for a broad range of network sizes 
TV. 

Previous works have investigated neural networks in the thermodynamic limit 
— 7- oo [11, 13], scaling the synaptic amplitudes to zero in order to arrive at 
analytical results. Equation (8) provides an alternative criterion to scale the synapses 
while keeping the dynamics comparable: The two terms in equation (8) depend 
differently on the feedback L. In order to maintain their ratio, we need to keep 
the population feedback L constant. The synaptic amplitude J (approximately 
oc w (A. 3)) hence needs to scale as J oc 1/A^. In addition, the response kernel 
h of each single neuron must remain unchanged, requiring the same working point, 
characterized by the mean /ij and fluctuations (Xj in the input into each cell. Constant 
mean directly follows from L = const., but the variance due to local input from other 
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neurons in the network decreases as 1/A^. To compensate, we supply each neuron 
with an additional external uncorrelated balanced noise whose variance appropriately 
increases with (as described in detail in Appendix E). Figure 3b shows that the 
shape of the covariance functions is invariant over a large range of network sizes A^, 
in particular the apparent time lag of inhibition behind excitation observed as the 
asymmetry in Figure 3b does not vanish in the limit N ^ oo. The magnitude of the 
covariance decreases as 1/A^ as expected from equation (8), because Kw = const. 
Previous studies have applied a different scaling J oc 1/ \/N in order to keep the 
variance at the input to each neuron constant [11]. Such a scaling increases the 
feedback on the network level L oc y/N and therefore changes the collective network 
state. 

Global properties of the network dynamics can be inferred by considering the 
spectrum of equation (7), those complex frequencies Zk at which the expression has 
a pole due to the function U{z). These poles are resonant modes of the network, 
where the real part ^{zk) denotes the damping of the mode, and the imaginary part 
'^{zk) is the oscillation frequency. A pole appears whenever Zk is a single root of 
U~^{zk) = H~^{—izk) — L = 0. With the Fourier representation H{uj) = ^jj^^ of 
the response kernel (5), the poles correspond to the spectrum of the delay differential 
equation Te^(t) = —y{t) + Ly{t — d), (cf. [29]) which describes the evolution of the 
population averaged activity. As shown in Appendix C, the location of the poles 
can be expressed by the branches k of the Lambert-iy function, the solution of 
WkC^^ = X [30], as 

^fc= -- + Vfc(^-e-) keNo. (9) 

The spectrum only depends on the population feedback L, the delay d and the 
effective time constant Tg of the neural response kernel. This explains why keeping 
L constant while scaling the network in Figure 3 yields shape invariant covariance 
functions. A typical spectrum is shown in Figure 4a as an inset, where each dot marks 
one of the solutions zj. of equation (9). The two principal branches of Wk are the 
modes with the largest real part ^{z^), and hence with the least damping, dominating 
the network dynamics. The remaining branches appear as conjugate pairs and their 
real parts are more negative, corresponding to stronger damping. Investigating the 
location of the principal branches therefore enables us to classify the dynamics in the 
network. Their dependence on the delay is shown in Figure 4a as a parametric plot in 
d. The point at which the two real principal solutions turn into a complex conjugate 
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Figure 4. Phase diagram determined by spectral analysis, (a) The two rightmost 
poles change with delay d (colors correspond to delays throughout). At d = 
0.753 ms (gray cross) (10) the poles become a conjugate pair, at c? = 6.88 ms 
(black crosses) (12) both poles have a zero real part, causing oscillations (Hopf 
bifurcation). Each dot in the inset represents a pole Zk (9) for delay d ~ 1 ms 
(two rightmost poles appear as one point), (b) Right: Phase diagram spanned 
by Te/d and feedback L. Onset of oscillations below the black curve (12) (Hopf 
bifurcation, black crosses in a), damped oscillations below the gray curve (10) (gray 
cross in a). Left: Oscillation frequency (12) at the Hopf bifurcation, (c) Average 
cross covariancc between excitatory neurons and theory (15) (black), (d) Average 
autocovariancc of excitatory neurons (5-pcak not shown). 



pair marks a transition from purely exponentially decaying dynamics to damped 
oscillations. This happens at sufficiently strong negative coupling L or sufficiently 
long delay d, precisely when the argument of Wk is smaller than — e^^ [30], leading 
to the condition 

L< -^e-^-\ (10) 

The gray cross marks this point in Figure 4a, the gray curve shows the corresponding 
relation of feedback and delay in the phase diagram Figure 4b. In the region 
below the curve, the dominant mode of fluctuations in the network is thus damped 
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oscillatory, whereas above the curve fluctuations are relaxing exponentially in time. 

For sufliciently long delay d the principal poles may assume positive real values, 
leading to ongoing oscillations, a Hopf bifurcation. The condition under which 
this happens can be derived from ^/"^(wcrit.) — -^^ = 0, as detailed in Appendix C. 
Equating the absolute values on both sides leads to the condition WcritTe = V L"^ — I: 
oscillations can only be sustained, if the negative population feedback is sufficiently 
strong L < —1. The oscillation frequency increases the stronger the negative 
feedback. The condition for the phases leads to the critical delay required for the 
onset of oscillations (see Appendix C for details) 

c^crit — arctanfV-^^ — 1) / s 

= ^ ^. (11) 

This relation is shown as the black curve in the phase diagram Figure 4b. The 
oscillatory frequency on the bifurcation line, at the onset of oscillations can be 
expressed as 

27r/crit.'i = oJcvit.d = 71 — arctan(A/ — 1), (12) 

which is shown in the left sub-panel of the phase diagram Figure 4b. Consequently, 
the oscillation frequency /crit. at the onset is between (4(icrit.)~^ and (2(icrit.)~^ , 
depending on the strength of the feedback L, with increasing negative feedback 
approaching /cnt. = (4(icrit.)~^ at the onset. 

Changing the synaptic delay homogeneously for all synapses in the network 
allows to observe the transition of the network from exponentially damped, to 
oscillatory damped, and finally to oscillatory dynamics. For a short delay of 
d = 0.5 ms the dynamics is dominated by the single real pole near — l/rg (brown 
dot in Figure 4a) and the covariance function is exponentially decaying (Figure 4c). 
Increasing the delay to d = 1 ms the principal poles split into a complex conjugate 
pair as the delay crosses the gray curve in Figure 4b so that side troughs become 
visible in the covariance function in Figure 4c. Further increasing the delay, the 
network approaches the point of oscillatory instability, where a Hopf bifurcation 
occurs, marked by black crosses in Figure 4a and the black curve in Figure 4b. 
The damping of oscillations decreases as the system approaches the bifurcation 
(Figure 4c). The auto covariance function of single spike trains, however, stays 
mostly unchanged (Figure 4d). This state is known as synchronous irregular activity 
[9], where the activity of a single neuron is irregular, but collectively the neurons 
participate in a global oscillation. 
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Figure 5. Composition of covariance functions, (a) Network echo caused by 
a spike sent at time t = sets in after one synaptic delay (d = 3 ms): red 
curve in c, inhibitory spike, blue curve in c, excitatory spike, (b) Correlated 
inputs from the network to a pair of neurons (black circles) cause covariance Cff 
between their outputs (green curve in c). (d)-(f) Covariance functions averaged 
over pairs of neurons (black dots, d: excitatory pairs, e: excitatory-inhibitory 
pairs, f : inhibitory pairs) and theory (15) (gray underlying curve). Inset shows the 
two components from c that add up. 
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If the network is not in the oscillatory state, all modes are damped in time, so 
all poles Zk of the function U{z) appearing in equation (7) lie in the left complex 
half plane, ^{zk) < 0; the function |f/(z)p has poles in both half-planes. We 
perform the Fourier transform to time domain using the residue theorem u{t) = 
U{z) e^* dz = YliZk&Et R'es(f/, Zi) e^*-*, where the integration path Et proceeds 
along the imaginary axis from —zoo to zoo and is then closed in infinity in the left 
half plane (for t > d) oi the right half plane (for < t < d) to ensure convergence, 
resulting in (see Appendix D for the detailed calculation) 



nit) = ^ \ 9(t - d)e^^('-'\ (13) 

^ (1 + ZkTe) d + Te 

The back transform of V{u) '= \U{u)\'^ proceeds along similar lines and results in 

""^^^ " ^ (l + ZfcTjd + Te (l-^fere)-Le-^'^" ^^^^ 

The population averaged covariance functions in the time domain then follows from 
equation (7) for t > as 



, , Kw ( 1 —q \ 

= (i -l ) 

(15) 

which is the central result of our work. Figure 5 shows the comparison of the theory 
(15) with the covariance functions obtained by direct simulation. The analytical 
expression unveils that the covariance functions are composed of two components: 
the first line in equation (15) has heterogeneous matrix elements and hence depends 
on the neuron types under consideration. Its origin is illustrated in Figure 5a: If 
one of the neurons emits a spike, as indicated, this impulse travels along its axon 
and reaches the target neurons after one synaptic delay d. Depending on the type 
of the source neuron, the impulse excites (synaptic amplitude J) or inhibits {—gJ) 
its targets. Its effect is therefore visible in the pairwise covariance function as a 
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positive (blue) or negative (red) deflection, respectively in Figure 5c. This deflection 
not only contains the direct synaptic effect, but also inflnitely many reverberations 
of the network, seen formally in equation (13). This expression is not proportional 
to the kernel h(t) directly, but is rather a series including the whole spectrum of the 
network. The shape of the spike echo consequently shows onsets of reverberations at 
integer multiples of the synaptic delay (Figure 5c), being transmitted over multiple 
synaptic connections. The contribution of the second line in equation (15) follows the 
intuitive argument illustrated in Figure 5b. The incoming activity from the network 
to a pair of neurons is correlated. As the input statistics is the same for each neuron, 
this contribution is identical for any pair of neurons (green curve in Figure 5c). The 
sum of both components results in the covariance functions shown in Figure 5d- 
f. The same analytical solution is shown for different delays in Figure 4c showing 
good agreement with direct simulation. For different sizes of simulated networks 
in Figure 3b-d the analytical expression (15) explains why the spike echo does not 
become negligible in the thermodynamic limit — )■ oo: for flxed population feedback 
L, both contributions in equation (15) scale as 1/A^, so the relative contribution of 
the echo stays the same. This also explains the apparent paradox (see Figure 1), that 
covariance functions in recurrent networks not only depend on the input statistics, 
but in addition the spike feedback causes a reverberating echo. The power spectrum 
of population-averages is dominated by pairwise covariances, explaining the different 
spectra observed in the excitatory and inhibitory population activity [15]. Scaling 
the network such as to keep the marginal statistics of single neurons constant, 
J oc w oc I/Vn [11, 13] changes the spectrum (9), because the feedback increases as 
L oc Vn which can ultimately lead to oscillations as shown in Figure 4c. 

3. Discussion 

The present work qualitatively explains certain features of the correlation structure 
of simultaneously recorded synaptic currents of two cells in vivo. Novel experimental 
techniques are able to separate contributions of excitatory and inhibitory inputs [20]. 
We calculate such covariances in a random network and show that the covariance 
between synaptic impulses decomposes into a linear combination of the covariance of 
the spiking activity and the autocovariance functions (see caption of Figure 6). Each 
synaptic impulse has a certain time course, here modeled as a flrst order low pass 
with time constant = 2 ms (1). The covariance between these filtered currents is 
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t (ms) 

Figure 6. Covariance between synaptic currents of a pair of neurons in a recurrent 
random network in analogy to in vivo experiments [20] . Covariance of excitatory 
contributions (blue, cj^j^ = q * {J^pKa^ + J^K^Ccc)), analogously between 
inhibitory contributions (red), between excitatory and inhibitory contribution 
(green, c/^/j = —q * (gJ^jK^Cei)), and Cj.j^ analogously (brown). Currents are 
filtered by an exponential kernel of time constant = 2 ms (1), leading to the 
filtering of the covariances by q{t) = ^e^'*'/'^'. 

shown in Figure 6, their temporal structure resembles those measured in cortex in 
vivo [20, Figure le,f]: covariances between afferents of the same type are monophasic 
and positive, while the covariances between excitatory and inhibitory afferents are 
biphasic and mostly negative. The lag reported between inhibitory and excitatory 
activity [20, Figure 2b], which was also observed in binary random networks [11, 16], 
is explained by the echo of the spike contributing to the covariance function. In 
contrast to previous work, we take the delayed and pulsed synaptic interaction into 
account. Without delays and with binary neurons [11, 16] the echo appears to be a 
time lag of inhibition with respect to excitation. 

Measuring membrane potential fluctuations during quiet wakefulness in the 
barrel cortex of mice [31], showed that correlations between inhibitory neurons are 
typically narrower than those between two excitatory neurons [31, their Figure 4A, 5B 
and Figure 5C,E]. These results qualitatively agree with our theory for covariances 
between the spiking activity, because fluctuations of the membrane potential are 
uniformly transferred to fluctuations of the instantaneous firing intensity. The direct 
measures of spiking activity [31, their Figure 6] confirm the asymmetric correlation 
between excitation and inhibition. The low correlation between excitatory neurons 
reported in that study may partly be due to the unnormalized, firing rate dependent 
measure and the low rates of excitatory neurons. 
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Oscillations in the 7 range (25 — 100 Hz) are ubiquitous in population measures 
of neural activity in humans, and have earlier been explained in networks of leaky 
integrate-and-fire model neurons [9] by the Hopf bifurcation induced by delayed 
negative feedback. For the regime of high noise here we uncover a simpler analytical 
condition for the onset (11) and frequency (12) of fast global oscillations. For lower 
noise, deviations of the non-linear leaky integrate-and-fire dynamics from the linear 
theory presented here are expected. 

From a physics viewpoint, neuronal networks unite several interesting properties. 
They do not reach thermodynamic equilibrium even in the stationary state, as no 
detailed balance condition holds between pairs of states of the system. This can 
be appreciated directly from the fact that a constant firing rate in the network 
amounts to a continuous flux of the neurons' state variables (voltages) over the 
threshold. Moreover, the interaction between pairs of neurons is directed, delayed, 
pulsed, and depends on the flux of the sending neuron's state variable at threshold. 
In contrast, pairwise interactions frequently studied in physics, like the Coulomb 
interaction or exchange interaction, can be expressed by a pair potential and are thus 
symmetric (undirected), instantaneous, and depend directly on the state variables 
(e.g. spatial coordinates or spins) of the pair of interacting particles. Non-equilibrium 
systems are at the heart of ubiquitous transport phenomena, like heat or electric 
conduction. Understanding fluctuations in such a system marks the starting point 
to infer macroscopic properties by the assertion of the fluctuation-dissipation theorem 
that connects microscopic fluctuations to macroscopic transport properties. Despite 
the non-equilibrium dynamics and the non- conservative pairwise interaction, in 
this manuscript we develop a simple analytical framework merely based on linear 
perturbation theory that explains time dependent covariance functions of the activity 
of pairs of integrate-and-fire model neurons in a recurrent random network. Formally 
our approach resembles the step from the kinetic Ising model near equilibrium to its 
non-equilibrium counterpart, the network of binary neurons [16]. A difference is 
the spiking interaction considered in our work, which led us to the describe each 
neuron in terms of the flux over threshold (spike train) rather than by its state 
variables (membrane voltage and synaptic current). In this respect, we follow the 
established mean- field approach for spiking neuronal systems [9, 14]. While this 
mean field approach, however, proceeds by assuming vanishing correlations to obtain 
the dynamics of the ensemble averaged activity, here we derive and solve the self- 
consistency equation for the pairwise averaged covariances of the microscopic system. 
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The typical time scale of covariance functions found here coincides with the time 
window of biological synaptic plasticity rules [21], so that non-trivial interactions of 
dynamics and structure are expected. It is our hope that the novel capability to 
resolve the temporal structure of covariances in spiking networks presented here 
proofs useful as a formal framework to further advance the theory of these correlated 
non-equilibrium systems and in particular serves as a further stepping stone in 
the endeavor to understand learning from the interplay of neuronal dynamics and 
synaptic plasticity. 
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Appendix A. Response kernel of the LIF model 

The response kernel kernel hij needs to be related to the dynamics of the neuron 
model (1). Here we present an approximation of this kernel which is sufficiently 
accurate to allow quantitative predictions, but yet simple enough to enable an 
analytical solution for the correlation structure. If the synaptic time constant is 
short Ts <^ Tm, the synaptic amplitude J can be thought of as the amplitude of 
the jump in the membrane potential V caused upon arrival of an incoming impulse. 
If correlations between incoming spike trains are sufficiently small, the ffist and 
second moments of the summed impulses TmYlj Jij^ji^ ~ ^) jii = TmJ2j 
and af = r^J^j --^ifj^ respectively, if the inputs' statistics can be approximated by 
Poisson processes of rate each. For small J and a high total rate, the system of 
differential equations (1) is hence approximated by a stochastic differential equation 
driven by a unit variance Gaussian white noise ^ 

dVi T / \ 

^■^^ ^ - Vi + L,i[t> 
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The stationary firing rate in this hmit is given by [24] 

r-r' =rr + r^V^ {F{ye) - F{yr)) (A.l) 
f{y) = ey\l + eriiv)) F (y) = J' f (y) dy 

with,,, = ^^^ + f./^ a = V2\ab\, 

with Riemann's zeta function (. The rate rj is the density of action potentials per 
time. The response of the firing density of the neuron i at time point t with respect 
to a point-like deflection of the afferent input Sj at time point t' defines the response 
kernel as the functional derivative 

\ lim - {si{t, {s(r) + eS{T - t')e,|r < t}) 



- s,(t,{s(r)|r<t}))^^^^ 
hij{t-t') =Wij h{t-t'). 

Here we used the homogeneity, namely the identical input statistics of each neuron 
2, leading to the same temporal shape h{t) independent of i and the stationarity, so 
that the kernel only depends on the time difference t — t'. We chose h{t) to have unit 
integral and defined Wij as the integral of the kernel. We determine the temporal 
integral of the kernel as 

/oo 
hij{t) dt (A. 2) 

-oo 

dr- 

The second equality holds because the integral of the impulse response equals the 
step response [32]. Further, a step in the density Sj corresponds to a step of r^-. Up to 
linear approximation the effect of the step in the rate r^ on the rate r^ can be expressed 
by the derivative considering the perturbation of the mean /ij and the variance al 
upon change of r^. Using equation (A.l) we note that by chain rule |^ = —'rl^-^- 

The latter derivative follows as = -^^ir^ + -r^ir^- The first derivative in both 

OTj aye arj clyr arj 

terms yields = y/TTTmf{yA) with yA € {y0,yr}- The second derivative evaluates 
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with yA = ^ + 



^ fl£ to = 




= -r, 



m 



inn + So 



together we obtain 



with a 



and /3 



yvr(rmri)2— (/(yg) - /(?/,.)) 
A/vr(r^r,)^— ( /(t/g) Jij-^^ 




fiVr) Ji 



2af 



) 



(A.3) 



Appendix B. Cross spectral matrix in frequency domain 

The autocovariance Aiu) in frequency domain {F{uj) = J^[f]{uj) = f{t)e~'^'^^ dt) 
has two different terms. The first term is a constant r due to the spiking with 
rate r resulting from the delta peak r6{t) in time domain. The second term is the 
continuous function adt), for example due to refractoriness of the neuron. Further 
follows from a{t) = a{—t) that A{u!) = A{—u!). For r > the covariance matrix 
fulfills the linear convolution equation (6). As this equation only holds for the 
positive half of the time axis, we cannot just apply the Fourier transform to obtain 
the solution. For negative time lags r < the covariance matrix is determined 
by the symmetry c(r) = c^(— r). Here we closely follow [17] and employ Wiener- 
Hopf theory [26] to derive an equation for the cross spectral matrix in frequency 
domain that has the desired symmetry and solves (6) simultaneously. To this end 
we introduce the auxiliary matrix b(r) = {h * (Mc)) (r) + Q(/i * a)(r) — c(r) for 
— oo < r < oo. Obviously, b(r) = for r > 0. Since the defining equation 
for b holds on the whole time axis, we may apply the Fourier transform to obtain 



We observe that QM"^ = MQ-^ is symmetric and with A{u) = A{—uj) the term 
proportional to \H{u)\'^ cancels on both sides, remaining with 



B(w) = H{u) (MC{u) + QA{u)) - C{u). Solving for C 

C{u) = (1 - MH{uj))'^ {H{uj)A{u)Q - B{u)) 
and using the symmetry C(ci;) = C'^{—u) we obtain the equation 



(B.l) 



{H{uj)A{uj)Q - B{uj)) (1 - M^H{~u)) 
(1 - MH{uj)) {H{-uj)A{-io)Q^ - B^(-w 



))• 



H{u)A{uj)Q + (1 - MH{uj)) B^{-u) 
H{-uj)A{uj)Q^ + B{uj) (1 - M^H{-uj)) . 



(B.2) 
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We next introduce D{uj) = H{u)A{u) which we spht into D{uj) = D^{uj) + D„{u), 
chosen such that d^{t) (in time domain) vanishes for t < and d^{t) vanishes 
for t > 0. Consequently the Fourier transforms of both terms may have poles in 
distinct complex half-planes: D^{uj) may only have poles in the upper half plane 
> and the function vanishes for lim|^|_j,oo,'3(u;)<o = 0, following from 
the definition of the Fourier integral. For D^{uj) the half planes are reversed. The 
analytical properties of Hiu) are thus similar to those of D^[bj), those of B[u:) are 
similar to D_{u:). We sort the terms in (B.2) such that the left hand side only 
contains terms that vanish at infinity in the lower half plane < 0, the right 

hand side those that vanish in infinity in the upper half plane '^{uj) > 

D+{uj)Q - D4-uj)Q^ + (1 - MH{uj)) B^(-a;) (B.3) 

= D+{-io)Q^ - D_(u;)Q + B{uj) (l - M^H{-uj)) . 

The left hand side consequently is analytic for ^^(w) < the right hand side is analytic 
for > 0, so (B.3) defines a function that is analytic on the whole complex plane 
and that vanishes at the border for — ?■ oo. Hence by Liouville's theorem it is 
and we can solve the right hand side of (B.3) for B 

B(a;) = (D_(cj)Q - D+{-uj)Qr^) (l - M^H{-uj)y^ . 
Inserted into (B.l) this yields with the definition P(a;) = (1 — MiJ(a;))~^ 

C{u) = P{cu) {H{lo)A{u)Q - {D_{cu)Q - D+(-a;)Q^) P^(-a;)) 
= P{uj) {{D+{uj) + D^u)) Q(l - M^H{-u)) 

-D^{uj)Q + D+{-u)Q^) P^{-uj) 
= P{uj) {D+{uj)Q + D+(-c^)Q^ - A{u)\H{uj)\^ClM^) P^{-uj). 

The latter expression can be brought to the form 

C{u) = D+{u)P{uj)Q + D+{-uj)QV{-u) 

+ {D+{uj)H{-uj) + D+{-uj)H{uj) - A{uj)\H{uj)\'^) x (B.4) 
X P(a;)MQ^P^(-u;), 

which has the advantage that the first two terms have poles in distinct half planes 
'^{uj) > 0, and '^{oo) < 0, respectively. This means these terms only contribute for 
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positive and negative times, respectively, the last term contributes for positive and 
negative times. 



Appendix C. Spectrum of the propagator 

— iud 

With the Fourier representation H{u) = of the delayed exponential kernel (5) 

the averaged cross spectrum (7) contains the two functions U{z) and V{z) = \U{z)\'^ 
defined on the complex frequency plane z = iu. These functions may exhibit 
poles. The function U in has a pole Zk whenever the denominator has a single 
root H~^{—iz) — L = which amounts to the condition (1 + ZkTe)e^'''^ = L. These 
complex frequencies can be expressed by the Lambert W function, the solution of 
We^ = X [30], by 



^1 + ZkTe)e 

[ h Zkdje"^ 



(C.l) 



+Zkd 



L — 



as 



Zk= + -Wk{L—e-e), 

Te d Te 

leading to (9). The Lambert Wk{x) function has infinitely many branches k [30]. 
The principal branch has two real solutions, if x > — e~^. The remaining branches 
appear in conjugate pairs. For x < —e~^ the principal solutions turn into a complex 
conjugate pair. This happens at sufficiently strong negative coupling L or long delays 
d 

d A. _i 
L — e^e < — e 

Te 

L < -e 

d 

The principal poles may assume positive real values, leading to oscillations. The 
condition under which this happens can be derived from (C.l). At the point of 
transition the pole can be written as 2 = iuJcxit.] it is a solution to (l+iaJcritTe)^*'^"'* '^ = 
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L. In order for this equation to be fulfilled, the absolute value and the phase 
must be identical on both sides. The equation for the absolute value requires 
1 + (wcrit.Te)^ = ■ This meaus there are only oscillatory solutions, if the magnitude 
of the feedback exceeds unity L < —1. Since the poles come in conjugate pairs, we 
can assume w.l.o.g. that Wcrit. > 0. The condition for the absolute value hence reads 



'^crit.'^'e 



VZ^^. (C.2) 



This is the frequency of oscillation at the onset of the Hopf bifurcation. For strong 
feedback |L| ^ 1 the frequency increases linearly with the magnitude of the feedback. 
The condition for the agreement of the phases reads Z(l + zWcritTe) + i^crit.c^ = 0, 
^° ^(i+Zcri! Te) = ~ ^[e'"cri!>] = tanwcrit.c?, which leads to tanwcrit.c? = -t^crit.Te. This 
equation has a solution in ^ < Ucvitd < n. In the limit of vanishing delay (i — ?■ the 
frequency goes to infinity, as the solution converges to oJcrit.d = |. This corresponds 
to the frequency fcnt. = j^- Inserting (C.2) leads to tan{^\^L'^ — 1) = — V — 1, 
which can be solved for the critical delay 



d 71 — arctan(v^L2 — 1) 

where we took care that the argument of the tangent is in [|, tt]. So with (C.2) and 
(C.3) the oscillatory frequency at the transition can be related to the synaptic delay 

as 



27ifcrit.d = ojd = 71 — a.Tctan{V L"^ — 1). 
Appendix D. Back transform by residue theorem 

In the non-oscillatory state all poles Zk (9) have a negative real part. The function 
U{z) = ((1 + zTf.)e^'^ — L)~^ in (7) then has all poles in the left complex half plane, 
3fJ(2fc) < V/c. We perform the Fourier back transform 

u{t) = — / U{iuj)e''^'duj (D.l) 
= — / U{z) e"* dz, 



E, 



replacing the integration path by a closed contour Et following the imaginary axis 
from —ioo to ioc. In order to ensure convergence of the integral, ioi t < d we need 
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^{z) > 0, so we close -E't<d in infinity within the right half-plane, where the integrand 
vanishes. Since there are no poles in the right half-plane, for t < d the path Et<:d 
does not enclose any poles, so u{t) = 0. For t > d the path Et>d must be closed in 
the left half-plane to ensure convergence of (D.l), so the residue theorem yields 

u{t) = e{t -d) Res{U, Zk) e"^*. (D.2) 

The residue can be calculated by linearizing the denominator of U {z^ + z) around z^ 



(1 + {zk + ^)re)e(^'=+")'^ - L = (1 + ^fer,)e"'=V'^ - L + ^r,e"'= V'^ 

= L{e"^ - 1) + zT^e^'^'^il + zd) + 0{z^) 
= z{Ld + Tee'^'^) + 0{z^), 

which yields 

Res(t/, Zk) = lim zU (zk + z) 
= lim 



2->0 z{Ld + Tet'^kd) Ld + TeC'^kd 
g— Zfcd 



(1 Zkre)d + Te 

where we used (C.l) in the last step. The poles of V{z) = U{z)U{—z) are located 
in both half-planes, consequently v(t) is nonzero on the whole time axis. Here we 
only calculate v{t) for positive times t > 0, because it follows for negative times by 
symmetry v{—t) = v{t). The path has to be closed in the left half-plane, where the 
poles Zk have the residues 

So applying (D.2) the functions u and v are 

oo ^ 

fe=0 
oo 



,{t) = ^Res(l^,2fc) (e(t)e"'=* + e(-t)e-"^*) 

A:=0 

oo . 

= V - t 

^ (1 + ZkTe) d + Te (1 - ^fcT^) - Lc'^'^ 
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The amplitude of the modes decrease with k. For all figures in the manuscript we 
truncated the series after k = 30. 

Appendix E. Simulation parameters used for figures 

All network simulations were performed using NEST [33]. The parameters of the 
leaky integrate-and-fire neuron model (1) throughout this work are Tm = 20 ms, 
Ts = 2 ms, Tr = 2 ms, Vg = 15 mV, Vr = mV. All simulations are performed 
with precise spike timing and time stepping of 0.1 ms [34]. Figure 1 and Figure 3- 
Figure 6 of the main text all consider recurrent random networks of excitatory and 
7 A inhibitory leaky integrate-and-fire model neurons receiving input from randomly 
drawn neurons in the network and external excitatory and inhibitory Poisson input, 
so that the first and second moments are /ij = 15 mV and erf = 10 mV, respectively. 
Unless stated explicitly, we use A(l + 7) = 10000 neurons except in Figure 3 where 
the number of neurons is given in the legend. Each neuron has K = pN incoming 
excitatory synapses with synaptic amplitude J independently and randomly drawn 
from the pool of excitatory neurons, and 7A = jpN inhibitory synapses with 
amplitude —gJ (homogeneous Erdos-Renyi random network with fixed in-degree), 
realizing a connection probability p = 0.1. Cross covariance functions are measured 
throughout as the covariance between two disjoint populations of 1000 neurons each 
taken from the indicated populations in the network. Correlation functions are 
evaluated with a time resolution of 0.1 ms. 

In Figure 3 we keep the feedback of the population rate constant L = Kwll — 
7(7) = const. Increasing the size of the network A the synaptic amplitude J (which 
is proportional to w in linear approximation) needs to scale as J = JqNq/N, where 
we chose Jq = 0.1 mV and Aq = 10000 here. The variance caused by local input 
from the network then decreases with increasing network size oc 1/A, while the local 
mean is constant because L = const. Each cell receives in addition uncorrelated 
external balanced Poisson input, adjusted to keep the mean /ij = 15 mV, and 
fluctuations cr^ = 10 mV constant. This is achieved by choosing the rates of 
the external excitatory (re,ext., amplitude Jext. = 0.1 mV) and inhibitory (ri^xt., 
amplitude —gJcxt.) inputs as 

'^e.ext. = '^e.O + '^bal '^i.ext. = '^bal/fi' (E.l) 

2 2 t2 

Atj — Mloc. -, CTi ~ CTlnr. ~ ''"m'^cO-'nxt. 
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where /xioc. = TmrKJ{l — 7(7) and a'(^^ = TmrK-Pil + 7(7^) are the mean and 
variance due to local input from other neurons of the network firing with rate r. 
From L = const, follows that also Kw = const., so that (8) predicts a scaling of 
the magnitude of the covariance functions in proportion to Other network 

parameters are d = 3 ms and g = 5. The firing rate in the network is r = 23.6 Hz. 

In Figure 4 and Figure 5 we use = 10^, g = 6, J = 0.1 mV, and the delay d 
as described in the captions, the remaining parameters are as in Figure 3. 

In Figure 6 we use N = 10^, J = 0.1 mV, g = 5, and d = 2 ms and the 
remaining parameters as in Figure 3. We obtain the filtered synaptic currents by 
filtering the spike trains with an exponential filter of time constant = 2 ms. This 
results in an effective filter for the cross covariances of = I^e-l*l/^^ The different 
contributions shown are C/^/^ = q* {J'^pKac + J'^K'^Cee), Ci.j. = q* {J'^pKai + J^K'^Ca), 
ci^j. = —q * (gJ'^jK'^Cci), and Cij^ = —q * {gJ'^'yK'^Cic), where * denotes the 
convolution. 

For Figure 2, we simulate two populations of = 1000 neurons each. Each 
neuron receives independent background activity from Poisson processes and in 
addition input from a common Poisson process with rate = 25 Hz causing in 
population 1 a positive synaptic amplitude of J and for population 2 a negative 
synaptic amplitude J (J is given on the x-axis). The synaptic amplitude of the 
background inputs is Je = 0.1 mV for an excitatory impulse and Ji = —0.5 mV for 
an inhibitory impulse. The rates of the excitatory and inhibitory background inputs 
are chosen so that the first and second moments /ij = Tm{Jc'^c + Jj^i + J^c) = 15 mV 
and af = Trn{Jlrc + Jiri + J r1) = 10 mV are independent of J. The spikes produced 
by each population are triggered to the arrival of an impulse in the common input 
and averaged over a duration of 10 s to obtain the impulse response. 
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